#include "Mesh/Cell_Traits.h"
#include "tool/libcell_tools_kd_node.h"
#include "tools/intrbtree.h"
#include "tools_node.h"
#include<stdio.h>
//#include<Mesh/Mesh_Frame.h>
//#include<Mesh/lib_cell_Iterator.h>
#include<Mesh_IO/Mesh_IO.h>
//#include<Algorithm/Algorithm1.h>
//#include<Algorithm/Algorithm2.h>
#include<tool/libcell_tools_formats.h>
#include<tool/libcell_tools_view.h>
//#include<Algorithm/Harmonic.h>
//#include <Algorithm/Conformal.h>
#include<tool/libcell_tools_octree.h>
#include<tool/libcell_tools_algorithm.h>
#include<Mesh/libcell_iterator.h>
#include <Algorithm/qem.h>
#include<Algorithm/points_cloud_algorithm.h>
#include<libcell_macros.h>
#include<time.h>

#define _ReadOff_   LibCell_ReadOff_
#define _ReadCell_  LibCell_ReadCell_
#define _ReadPly_ LibCell_ReadPly_

#define _WriteCell_ LibCell_WriteCell_

void test_tools_algorithm()
{
    Mesh m,m1;
    Mesh_init(&m);Mesh_init(&m1);

    _ReadOff_(&m1,"face.off",3);
    _ReadCell_(&m,"cube.cell");
    printf("num f :%d\n",m.num_v(&m));
    template_v* v=m.get_vertexp(&m,0);
    Node* node=NULL;
    for(int i=0;i<v->cellsn.size;i++)
    {
        node=node_overlying(node,v->cellsn.extra_data[i]);
    }
    // node_copy(v->faces);
    //printf("node size:%d\n",node_size(node));
    Int_RB_Tree*tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree);

    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    Mesh_cut_along_the_curvef(&m,node,tree,tree1);

    printf("num f :%d\n",m.num_v(&m));
    for(auto vit=m.v_begin(&m);vit!=m.v_end(&m);vit++)
    {
        printf("%lf %lf %lf\n",quote(vit)->point[0],
            quote(vit)->point[1],
            quote(vit)->point[2]);
    }

    int_rb_tree_free_value(tree);
    int_rb_tree_free(tree);
    int_rb_tree_free(tree1);
    Mesh_free(&m); Mesh_free(&m1);
}
void test_read_ply()
{
    Mesh m;
    Mesh_init(&m);
    _ReadPly_(&m,"bunny.ply" );
    Mesh_free(&m);
}
void my_sscanf(char * str,int* i1,int* i2,int *i3)
{
    char str1[100]={NULL};
    int temp_i=0;
    sscanf(&str[temp_i],"%[^/]",str1);
    //printf("%s ",str1);
    sscanf(str1,"%d",i1);
    temp_i=temp_i+strlen(str1)+1;
    int i=sscanf(&str[temp_i],"%[^/]",str1);
    if(i==1)
    {
        temp_i=temp_i+strlen(str1)+1;
        //printf("%s %d ",str1,i);
        sscanf(str1,"%d",i2);
    }
    else
    {
        *i2=-1;
        temp_i=temp_i+1;
       // printf(" ");
    }
    i=sscanf(&str[temp_i],"%[^/]",str1);
    if(i==1)
    {
        sscanf(str1,"%d",i3);
        // printf("%s \n",str1);
    }
    else
    {
        *i3=-1;
       // printf("\n");
    }
}
double piecewise_map(double x);


double compute_distance(void*,double* );

void* fitting_plane(double**,int len);

void* compute_plane(double**,int len);




int * get_random_sequence(int len,int n)
{
    int* re=(int*)malloc(sizeof(int)*n);
    double initx=(rand()%900000)/900000.0;

    for(int i=0;i<n;i++)
    {
        //re[i]=(rand()%900000)/900000.0;
        initx= piecewise_map(initx);
        re[i]=floor( initx*(len-i));

    }

    for(int i=1;i<n;i++)
    {
        for(int j=0;j<i;j++)
        {
            if(re[i]>=re[j])
            {
                re[i]++;
            }
        }
    }


    return re;
}


double piecewise_map(double x)
{
    if(x<=0.00000001)
    {
        return 0.00001;
    }
    else if(x<0.30)
    {
        return x/0.30;
    }
    else if(x<0.5)
    {
        return (x-0.30)/(0.5-0.30);
    }
    else if(x<0.70)
    {
        return (1-0.30-x)/(0.5-0.30);
    }
    else{
        return (1-x)/0.30;
    }
    return 0;

}

void compute_ransac(double (*c)(void*,double*) ,void* (* f)(double**,int),
    void* (*cp)(double**,int),double**S,int len,int len1,double threshold)
{

    double** s= (double**)malloc(sizeof(double*)*len1);
    int*seq=get_random_sequence(len,len1);
    for(int i=0;i<len1;i++)
    {
        s[i]=S[seq[i]];
    }
    free(seq);


    void * p=f(s,len1);
    int sum=0;

    for(int i=0;i<len;i++)
    {
        if(c(p,S[i])<threshold)
        {
            sum++;
        }
    }


    free(p);

    free(s);
    // for()
    // {

    // }
}
#include<Math/LB_Math.h>
#include <Matrix/LB_Matrix.h>

static void gauss_bonnet(Mesh* m)
{
    double re=0,re1=0;
    LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,3,3);
    double* data=(double*)(mat->data);

    for(auto vit=m->v_begin(m);vit!=m->v_end(m);vit++)
    {
        double sum1=0;
        for(auto vfit=m->vf_begin(m,*quote(vit));vfit!=m->vf_end(m,*quote(vit));vfit++)
        {
            double data1[4][3]={0};
            template_c* c1=NULL,*c2=NULL;

            if(quote(vfit)->halffaces[0]->vertices[0]==quote(vit))
            {
                c1=quote(vfit)->halffaces[0]->cell;
                c2=quote(vfit)->halffaces[1]->cell;
            }
            else
            {
                c2=quote(vfit)->halffaces[0]->cell;
                c1=quote(vfit)->halffaces[1]->cell;
            }
            for(int i=0;i<c1->vertices_size;i++)
            {
                data1[i][0]=c1->vertices[i]->point[0];
                data1[i][1]=c1->vertices[i]->point[1];
                data1[i][2]=c1->vertices[i]->point[2];
                if(c1->vertices[i]!=quote(vit) )
                {
                    double p[3]={data1[i][0]-quote(vit)->point[0],
                        data1[i][1]-quote(vit)->point[1],data1[i][2]-quote(vit)->point[2] };
                    normalize(p,3);
                    data1[i][0]=quote(vit)->point[0]+p[0];
                    data1[i][1]=quote(vit)->point[1]+p[1];
                    data1[i][2]=quote(vit)->point[2]+p[2];
                }
            }
            for(int i=0;i<c2->vertices_size;i++)
            {
                if(m->vertex_is_in_cell(c2->vertices[i],c1))
                {
                    continue;
                }
                data1[3][0]=c2->vertices[i]->point[0];
                data1[3][1]=c2->vertices[i]->point[1];
                data1[3][2]=c2->vertices[i]->point[2];
                double p[3]={data1[3][0]-quote(vit)->point[0],
                        data1[3][1]-quote(vit)->point[1],data1[3][2]-quote(vit)->point[2] };
                normalize(p,3);
                data1[3][0]=quote(vit)->point[0]+p[0];
                data1[3][1]=quote(vit)->point[1]+p[1];
                data1[3][2]=quote(vit)->point[2]+p[2];
            }
            for(int i=0;i<3;i++)
            {
                for(int j=0;j<3;j++)
                {
                    data[i*3+j]=data1[i+1][j]-data1[0][j];
                }
            }

            double* det=(double*)(mat->det(mat));

            sum1+=(*det);

            SAFE_FREE(det);
        }
        printf("sum1 :%lf \n",sum1);

        re1+=sum1;

        double sum=0;
        for(auto vcit=m->vc_begin(m,*quote(vit));vcit!=m->vc_end(m,*quote(vit));vcit++)
        {
            template_v* vs[2]={NULL};
            int k=0;
            for(int i=0;i<quote(vcit)->vertices_size;i++)
            {
                if(quote(vcit)->vertices[i]!=quote(vit) )
                {
                    vs[k]=quote(vcit)->vertices[i];
                    k++;
                }
            }

            double dir1[3]={vs[0]->point[0]-quote(vit)->point[0],
                vs[0]->point[1]-quote(vit)->point[1], vs[0]->point[2]-quote(vit)->point[2] };
            double dir2[3]={vs[1]->point[0]-quote(vit)->point[0],
                vs[1]->point[1]-quote(vit)->point[1], vs[1]->point[2]-quote(vit)->point[2] };

            normalize(dir1,3);normalize(dir2,3);
            sum+=acos(inner_product(dir1,dir2,3));
        }
        re+=(2*M_PI-sum);
    }
    printf("gauss bonnet gongshi:%lf %lf\n",re,re1 ) ;

}

void test_read_obj()
{
    Mesh mesh;
    Mesh_init(&mesh);
    printf("begin read obj\n");

    _ReadObj_ (&mesh,"cat.obj");



    Mesh_free(&mesh);
}

void test_qem()
{
    printf("begin qem\n");
    Mesh mesh;
    Mesh_init(&mesh);
    _ReadOff_(&mesh,"face.off",3);
    printf("num v :%d num f:%d num c:%d\n",mesh.num_v(&mesh),mesh.num_f(&mesh),mesh.num_c(&mesh));


    Int_RB_Tree* tree=qem_init_qi(&mesh);

    CHeap*ch =compute_qem_all_edge_position_and_energy(&mesh,tree);
    for(int i=0;i<1000&&ch->size>0;i++)
    {
        //printf("onc ch size :%d\n",ch->size);
        void* value=pop_cheap(ch);

        Edge_Energy*ee=((Edge_Energy*)value);
        //printf("ee fid:%d size:%d\n",ee->fid,ch->size);
        if(qem_collapse_edge(&mesh,ee,tree,ch)==NULL)
        {
            i--;
        }
        //free(value);
        edge_energy_free(ee);
        //free(value);
    }

    printf("num v :%d num f:%d num c:%d\n",mesh.num_v(&mesh),mesh.num_f(&mesh),mesh.num_c(&mesh));
    _WriteCell_(&mesh,"qemmesh.cell");

    int_rb_tree_free(tree);
    Mesh_free(&mesh);
    printf("end qem\n");

}
// piont :-3.161449 -26.890598 14.075047
// piont :-3.299338 -27.011568 15.062319
// piont :-3.200117 -26.891392 12.861887
void test_triangle()
{
    double p1[3]={-3.161449, -26.890598, 14.075047};
    double p2[3]={-3.299338, -27.011568, 15.062319};
    double p3[3]={-3.200117, -26.891392, 12.861887};
    double dis= distance_of_two_points(p1,p3,3);
    double dis1=distance_of_two_points(p1,p2,3);
    double dis2=distance_of_two_points(p3,p2,3);
    double angle1 = acos((dis2*dis2 + dis1*dis1 - dis*dis) / (2 * dis1 * dis2));
    printf("angle1:%lf\n",angle1);
}
void test_external()
{
    printf("begin external\n");
    Mesh m,m1;
    Mesh_init(&m);Mesh_init(&m1);
    _ReadOff_(&m1,"bone.off",3);
    _ReadOff_(&m,"pmout.off",3);
    m.external_cells_init_(&m);
    for(Node* nit=m.external_cells;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c*c=(template_c*)(nit->value);
        if(c->vertices_size<60){continue;}
        int size=c->halffaces_size;
        for(int i=0;i<size;i++)
        {
            printf("%d %d\n",c->halffaces[i]->vertices[0]->id,
               c->halffaces[i]->vertices[1]->id);
            // pm->create_halfface(pm,c->halffaces[i]->face,
            //     c->halffaces[i]->vertices,c->halffaces[i]->vertices_size);
        }
        // pm->create_cellf(pm,c->halffaces,size);

        break;
    }
    // m.external_cells_free_(&m);
    _WriteOff_(&m,"pmout_test.off");

    // Mesh_free(&m1);Mesh_free(&m);
    printf("end read pmout off\n");
}

void my_test_fill_loop()
{
    Mesh mesh;
    Mesh_init(&mesh);
    _ReadOff_(&mesh,"右7牙颈线.off",3);

    int len=mesh.num_v(&mesh);
    template_v**vs=(template_v**)malloc(sizeof(template_v*)*len);
    double** vps=(double**)malloc(sizeof(double*)*len);
    memset(vps,0,sizeof(double*)*len);
    memset(vs,0,sizeof(template_v*)*len);
    int i=0;
    for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
    {
        vs[i]=quote(vit);
        vps[i]=quote(vit)->point;
        i++;
    }

    double normal[3]={0};
    compute_plane_normal_of_jordan_curve(vps,len,normal );
    printf("herer1\n");
    simple_fill_loop_hole_with_points(&mesh,
        vs,len,0.4);


    _WriteOff_(&mesh,"yajing.off");
    // Mesh_free(&mesh);
    printf("end yajing\n");

}
void test_external_cells()
{
    Mesh m;
    Mesh_init(&m);
    _ReadCell_(&m,"cube.cell");
    m.external_cells_init_(&m);
    for(Node* nit=m.external_cells;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c*c=(template_c*)(nit->value);
        for(int i=0;i<c->vertices_size;i++)
        {
            printf("%d ",c->vertices[i]->id);
        }
        printf("end once extern cell\n");
    }
}

// void test_fill_one_mesh(Mesh* m,Mesh*m1)
// {
//     template_v**vs=(template_v**)malloc(sizeof(template_v*)*m1->num_v(m1));
//     int i=0;
//     for(auto vit=m1->v_begin(m1);vit!=m1->v_end(m1);i++,vit++)
//     {
//         vs[i]=quote(vit);

//     }
//     int len=m->num_v(m);

//     simple_fill_loop_hole_with_points( m,
//         vs,len,0.3);

// }
//
Node* test_get_one_bound(Mesh* m,template_c*max_c,Mesh* ya)
{
    KD_Node* kdtree=create_kd_tree_from_mesh(ya);

    Int_RB_Tree* treev=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(treev);

    for(int i=0;i<max_c->vertices_size;i++)
    {
        treev->insert(treev,max_c->vertices[i]->id,max_c->vertices[i]);
    }
    Node* ren=NULL;
    for(int i=0;i<max_c->vertices_size;i++)
    {
        LB_Vector_voidptr vec= Mesh_vv_beginn(m,*(max_c->vertices[i]));
        template_v*v=NULL;
        double dis=0;
        for(int j=0;j<vec.size;j++)
        {
            template_v* v1=(template_v*)(*vec.at(&vec,j));
            if(treev->find(treev,v1->id)!=NULL)
            {
                continue;
            }
            double temp_dis=distance_of_two_points(v1->point,
                max_c->vertices[i]->point,3);
            if(v==NULL||temp_dis<dis){
                v=v1;
                dis=temp_dis;
            }
        }
        if(v==NULL){
            continue;
        }
        KD_Node*kdnode= kd_tree_find_nearest(v,kdtree);
        if(distance_of_two_points(v->point,kdnode->value->point,3)<0.2)
        {
            ren=node_overlying(ren, max_c->vertices[i]);
        }

        vec.clear(&vec);
    }


    int_rb_tree_free(treev);
    free_kdnode(kdtree);
    // template_v** re=NULL;
    // return re;
    return ren;
}
static void compute_cell_area_normal1(template_c*c,double* normal)
{
    if(c==NULL)
    {
        return ;
    }
    template_hf* hf=(template_hf*)(c->halffaces[0]);
    template_v* v=NULL;
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]!=hf->vertices[0]&&c->vertices[i]!=hf->vertices[1])
        {
            v=c->vertices[i];
            break;
        }
    }

    if(v==NULL)
    {return ;}

    double p1[3]={hf->vertices[0]->point[0]-v->point[0],hf->vertices[0]->point[1]-v->point[1],
        hf->vertices[0]->point[2]-v->point[2]};

    double p2[3]={hf->vertices[1]->point[0]-v->point[0],hf->vertices[1]->point[1]-v->point[1],
        hf->vertices[1]->point[2]-v->point[2]};
    out_productn(p1, p2, normal);
    // double*re=out_product(p1,p2);
    // normalize(normal,3);
    // return re;
}


static void compute_vertex_normal1(template_v*v,double*normal)
{
    // LB_Vector_voidptr vec=Mesh_vv_beginn(NULL,*v);
    memset(normal, 0, sizeof(double)*3);
    for(int i=0;i<v->cellsn.size;i++)
    {
        double normal1[3]={0};
        template_c* c=(template_c*)(*(lb_vector1_at_voidptr(&(v->cellsn),i)) );
        compute_cell_area_normal1(c,normal1);
        normal[0]+=normal1[0];
        normal[1]+=normal1[1];
        normal[2]+=normal1[2];
    }
    normalize(normal, 3);
}
void test_fill_one_hole(Mesh* m,template_c*max_c,Mesh* ya)
{
    Node* bound_node1= test_get_one_bound(m,max_c,ya);
    int len1=node_size(bound_node1);
    template_v** bound1=(template_v**)malloc(sizeof(template_v*)*len1);
    int i=0;
    double ** ps=(double**)malloc(sizeof(double*)*len1);
    for(Node* nit=node_reverse(bound_node1);nit!=NULL;
        i++,nit=(Node*)(nit->Prev))
    {
        template_v*v=(template_v*)(nit->value);
        bound1[i]=v;
        ps[i]=v->point;
    }
    free_node(bound_node1);
    double re[3]={0};

    compute_plane_normal_of_jordan_curve(ps,len1, re );
    free(ps);
    Node* nv_node=simple_fill_loop_hole_with_points(m,bound1,len1,0.3);

    KD_Node* kdtree=create_kd_tree_from_node(nv_node);

    points_cloud_compute_normalr(nv_node,kdtree,0.5);
    for(Node* nit=nv_node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v*v=(template_v*)(nit->value);
        double* normal=(double*)(v->prop);
        if(inner_product(normal,re,3)<0)
        {
            normal[0]=-normal[0];
            normal[1]=-normal[1];
            normal[2]=-normal[2];
        }
    }

    free_node(nv_node);
    free_kdnode(kdtree);

}

int _WriteNOff_(template_m* m,char const *filename)
{
    // printf("hrerebegin\n");
    if((quote(m->v_begin(m)))->point_size!=3)
    {
        printf("not 3\n");
        return 0;
    }
#ifdef _WIN32
    FILE *outfile;
    fopen_s(&outfile,filename,"w");
#else
    FILE *outfile=fopen(filename,"w");
#endif
    if(!outfile)
    {
        printf("cant open this file\n");
        return 0;
    }
    fprintf(outfile,"NOFF\n");
    fprintf(outfile,"%d %d %d\n",m->num_v(m),m->num_c(m),0);
    Int_RB_Tree* treev=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(treev);
    int *temp_v_id=(int*)malloc(sizeof(int)*m->num_v(m));
    int temp_id=0;
    for(auto v_it=m->v_begin(m);v_it!=m->v_end(m);v_it++)
    {
        temp_v_id[temp_id]=temp_id;
        treev->insert(treev,quote(v_it)->id,(&temp_v_id[temp_id]));
        double*normal=(double*)(quote(v_it)->prop);
        //((template_v*)(quote(v_it)))->prop=(void*)(&temp_v_id[temp_id]);
        for(int j=0;j<3;j++)
        {
            fprintf(outfile,"%.7lf ",((template_v*)(quote(v_it)))->point[j]);
        }
        if(normal==NULL)
        {
            fprintf(outfile,"%.7lf %.7lf %.7lf ",1.0,0.0,0.0);
        }
        else {
            fprintf(outfile,"%.7lf %.7lf %.7lf ",normal[0],normal[1],normal[2]);
        }
        fprintf(outfile,"\n");
        temp_id++;
    }
    int i=0;
    for(auto c_it=m->c_begin(m);c_it!=m->c_end(m);c_it++)
    {
        fprintf(outfile,"%d ",((template_c*)(quote(c_it)))->vertices_size);
        for(auto cv_it=m->cv_begin(m,*((template_c*)(quote(c_it))));cv_it!=m->cv_end(m,*((template_c*)(quote(c_it))));cv_it++)
        {
            int* vvluae=(int*)(treev->find(treev,(*cv_it).id));
            fprintf(outfile,"%d ",*vvluae);
        }
        fprintf(outfile,"\n");
        i++;
    }
    free(temp_v_id);
    // m->reset_v_prop(m,0);

 //   printf("herere libo ***********************************\n");

   int_rb_tree_free(treev);
   fclose(outfile);
   return 1;
}

void test_fill_hole()
{
    Mesh m;
    Mesh_init(&m);
    _ReadOff_(&m,"大连接体轮廓顶点/底片.off",3);
    for(auto vit=m.v_begin(&m);vit!=m.v_end(&m);vit++)
    {
        double * normal=(double*)malloc(sizeof(double)*3);
        compute_vertex_normal1(quote(vit),normal);
        quote(vit)->prop=normal;

    }


    m.external_cells_init_(&m);
    template_c* max_c=NULL;
    for(Node* nit=m.external_cells;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        if(max_c==NULL||c->vertices_size>max_c->vertices_size)
        {
            max_c=c;
        }
    }


    Mesh mp1;
    Mesh_init(&mp1);
    _ReadOff_(&mp1,"大连接体轮廓顶点/局部1牙齿.off",3);


    Mesh mp2;
    Mesh_init(&mp2);
    _ReadOff_(&mp2,"大连接体轮廓顶点/局部2牙齿.off",3);

    Mesh mp3;
    Mesh_init(&mp3);
    _ReadOff_(&mp3,"大连接体轮廓顶点/局部3牙齿.off",3);
    test_fill_one_hole(&m,max_c,&mp1);
    test_fill_one_hole(&m,max_c,&mp2);
    test_fill_one_hole(&m,max_c,&mp3);

    // Node* bound_node1= test_get_one_bound(&m,max_c,&mp1);
    // int len1=node_size(bound_node1);
    // template_v** bound1=(template_v**)malloc(sizeof(template_v*)*len1);
    // int i=0;
    // double ** ps=(double**)malloc(sizeof(double*)*len1);
    // for(Node* nit=node_reverse(bound_node1);nit!=NULL;
    //     i++,nit=(Node*)(nit->Prev))
    // {
    //     template_v*v=(template_v*)(nit->value);
    //     bound1[i]=v;
    //     ps[i]=v->point;
    // }

    // double re[3]={0};

    // compute_plane_normal_of_jordan_curve(ps,len1, re );

    // Node* nv_node=simple_fill_loop_hole_with_points(&m,bound1,len1,0.3);

    // KD_Node* kdtree=create_kd_tree_from_node(nv_node);

    // points_cloud_compute_normalr(nv_node,kdtree,0.5);


    // free_node(nv_node);
    // free_kdnode(kdtree);
    // _WriteOff_(&m,"outpoints.off");
    _WriteNOff_(&m,"outpoints.off");




    // void simple_fill_loop_hole_with_points(Mesh* m,
    //     template_v**vs,int len,0.3);
}

int main(int argc,char**argv)
{

    printf("sizeo f vector:%d\n",sizeof(LB_Vector1_voidptr));
    test_fill_hole();
    test_external_cells();

    printf("Size of LB_Vector1_voidptr: %zu bytes\n", sizeof(Cell));
    clock_t t=clock();

    my_test_fill_loop();
    test_external();
    printf("time si %.5f \n",(double)(clock()-t)/CLOCKS_PER_SEC );

    int inint=0;
    scanf("%d",&inint);
    return 0;
    test_triangle();
    test_qem();
    //return 1;

    double initx=0.121;

    for(int i=0;i<600;i++)
    {
        initx= piecewise_map(initx);
        printf("%lf ",initx);
    }
    test_read_obj();

    test_read_ply();
    test_tools_algorithm();
   //  FILE * infile=fopen("test.txt","r");
   //  int t=0;
   //  fscanf(infile,"%d",&t);
   //  printf("%d\n",t );
   //  char c;
   //  fscanf(infile,"%c",&c);
   //  printf("%c\n",c);
   //  fscanf(infile,"%c",&c);
   //  printf("%c\n",c);
   //  fscanf(infile,"%d",&t);
   //  printf("%d\n",t);
    printf("begin read cell\n");
    for(int i=0;i<1;i++)
    {
        Mesh m;
        Mesh_init(&m);
        _ReadCell_ (&m,"bone.cell");

        // _ReadOff_(&m,"bunny.off",3);
        // _WriteCell_(&m,"bunny.cell");
       Mesh_free(&m);
    }
    printf("end read cell\n");



    char s[200]="f 135//135 134/228 222/229/222 223/230/223";

    char str[100]={NULL};

    int temp_i=2;
    while(sscanf(&s[temp_i],"%s",str)!=EOF)
    {
        printf("%s %d\n",str,strlen(str));
        int ii1=0,ii2=0,ii3=0;
        char ss1[100]={NULL},ss2[100]={NULL};
        //sscanf(str,"%d[^/]%d[^/]",&ii1,&ii2);
        my_sscanf(str,&ii1,&ii2,&ii3);
        printf("%d %d %d\n",ii1,ii2,ii3);
        temp_i=temp_i+strlen(str)+1;
        memset(str,0,sizeof(char)*100);
    }

    // Mesh m;
    // Mesh_init(&m);
    // printf("%d\n",m.num_c(&m) );

    // objtooff("bone_scan1.obj");

 /*   mpf_set_default_prec(200);
    Mesh mesh,mesh1;
    Mesh_init(&mesh);Mesh_init(&mesh1);

    _ReadCell_(&mesh,"cube.cell");
    template_v v0=mesh.get_vertex(&mesh,0);
    template_f f0=mesh.get_face(&mesh,0);
    // mesh.delete_vertex(&mesh,v0,true);
    mesh.delete_face(&mesh,f0,true);
    double M[3][3],M1[3][2];
    M[0][0]=21.0;M[0][1]=-2.41;M[0][2]=-212.97;
    M[1][0]=1.001;M[1][1]=-0.51;M[1][2]=-2.4;
    M[2][0]=2.9;M[2][1]=92.25;M[2][2]=12.23;
    M1[0][0]=1;M1[0][1]=0;
    M1[1][0]=1;M1[1][1]=0;
    M1[2][0]=0;M1[2][1]=1;
    Tensors_Algebra_System*tas=(Tensors_Algebra_System*)malloc(sizeof(Tensors_Algebra_System));
    Tensors_Algebra_System_mpf_init(tas,3);
    Tensors_Algebra_System* tas1=(Tensors_Algebra_System*)malloc(sizeof(Tensors_Algebra_System));
    Tensors_Algebra_System_double_init(tas1,3);
    Tensor* t=tas1->T_create_matrix(tas1,(double*)M,3,3);
    //tas1->print_T(t);
    Tensor* t1=tas1->T_create_matrix(tas1,(double*)M1,3,2);
    tas1->print_T(t1);
    Tensor* t2=W_Tensor_Contraction(tas1,t,t1,1,0);
    tas1->print_T(t2);

    //tensor_double_print_self(t);

    //test_area();
    //test_convex();
    test_delauny();
    */
    //test_cell();
    //printf("%f\n",area_simplex_double(M,3,3));
    //Tensor*t=Anti_tensor_mpf_from_v(tas,M,3,3);
    //tensor_mpf_print_self(t);
    //gmp_printf("%.Ff\n",(__mpf_struct*)tas->T_norm(tas,t));

   // mesh.printself(&mesh);
 //   printf("v:%d,f:%d ,c:%d\n",mesh1.num_v(&mesh1),mesh1.num_f(&mesh1),mesh1.num_c(&mesh1));

   /* for(auto it=mesh.cells.begin();it!=mesh.cells.end();it++)
    {
        //printf("d");

    }*/
   // _ReadCell_(&mesh,"hand.cell");

    //objtooff("cat.obj");
    //pcdtooff("mypcd.pcd");
    //pcdtooff("tree1.pcd");
    // Mesh pm;
    // Mesh_init(&pm);
    // _ReadPly_(&pm, "luozhou_cloud.ply");
    //  Mesh pm1;
    // Mesh_init(&pm1);
    // _ReadObj_(&pm1, "checkerboard.obj");

    Mesh pm2;
    Mesh_init(&pm2);
    printf("begin read off\n");
    _ReadOff_(&pm2,"lawen.off",3);

    // gauss_bonnet(&pm2);

    //_ReadPly_
    //_ReadPly_(&pm, "bunny.ply");


    printf("end\n");
    return 0;
}
